knitr::opts_chunk$set(
    echo = TRUE,
    warning = FALSE,
    comment = "##",
    prompt = TRUE,
    tidy = TRUE,
    tidy.opts = list(width.cutoff = 75),
    fig.path = "img/"
)

Key Concepts

Geographic Coordinate Systems refer to data that is defined by a 3-D surface and measured in latitude and longitude. The longitude is the east or west location of a point, and the latitude is the north or south location of a point. Latitude and longitude coordinates thus define locations in two or three-dimensional space. An example of a Geographic Coordinate System would be “WGS84”. You may also encounter the word “Datum” when you work with GIS. Essentially a Datum provides a frame of reference for measureing locations on the surface of the earth, so “Datum” and “Geographic Coordinate System” are sometimes used interchangeably.

Projected Coordinate Systems refers to data that is defined by a flat 2-D surface and can be measured in units of meters and feet. The most important example is the The Universal Transverse Mercator projection (UTM). This is a cylindrical projection of the world’s layout on a map that touches at the equator. A projection is a mathematical transformation of the globe onto some other surface, in this case, a grid. UTM is the typical projection used in GIS and spatial analyses.

Shapefiles (.shp, .dbf, .shx, .prj, .sbn) store nontopological geometry and attribute information for the spatial features in a data set. The geometry for a feature is stored as a shape comprising a set of vector coordinates. Shapefiles can support point, line, and area features. One can create a shapefile on their map by downloading from an external source, digitizing, or programming. Shapefiles can be projected on a map in layers and manipulated.

Attributes are feature characteristics that can be plotted on a map.

Vector data models use discrete elements such as points, lines, and polygons to represent the geometry of real-world entities.

Point single point location, such as a GPS reading or a geocoded address

Line set of ordered points, connected by straight line segments

Polygon area, marked by one or more enclosing lines, possibly containing holes

Raster data models are the natural means to represent continuous spatial features or phenomenon. They define the world as a regular set of cells in a grid pattern that are square and evenly spaced in the x and y directions.

Grid collection of points or rectangular cells, organised in a regular lattice (raster data)

Before we take you on our expedition to learn how to analyze home range data, be sure to install the following packages in R: {maptools}, {sp}, {rgdal}, {zoom}, {adehabitatHR}, {curl}, {ggplot2}.

> library(maptools)
## Loading required package: sp
## Checking rgeos availability: TRUE
> library(sp)
> library(rgdal)
## rgdal: version: 1.2-4, (SVN revision 643)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.0.1, released 2015/09/15
##  Path to GDAL shared files: C:/Users/Silvy/Documents/R/win-library/3.3/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.2, 08 September 2015, [PJ_VERSION: 492]
##  Path to PROJ.4 shared files: C:/Users/Silvy/Documents/R/win-library/3.3/rgdal/proj
##  Linking to sp version: 1.2-3
> library(adehabitatHR)
## Loading required package: deldir
## deldir 0.1-12
## Loading required package: ade4
## Loading required package: adehabitatMA
## Loading required package: adehabitatLT
## Loading required package: CircStats
## Loading required package: MASS
## Loading required package: boot
> library(ggplot2)
> library(curl)

Section 2: Whatever the title was…

Part 1: Importing data

Importing and projecting a Shapefile

For this next section we will be looking at data from three groups of titi monkeys (Callicebus discolor) that are studied at the Tiputini Biodiversity Station in Ecuador.
Callicebus discolor

These data have been collected with a GPS. Every twenty seconds the GPS collects a ranging point automatically, which consists of the x and y coordinates for the location, altitude of the location, and a timestamp. From that, we (read: Tony) have created a .csv file that contains a subset of these points, twenty minutes apart. For today’s purposes we will use a subset of data from a one-year period: July 2014 to June 2015.

Let’s first load in a polygon of the trail sytem at research station. This can be done with the readShapeLines() command in the maptools package.

> tbs <- readShapeLines("C:/Users/Silvy/Documents/Austin/Tiputini/GIS/trails.shp")
> # tbs <- readShapeLines('~/Desktop/GIS+R_Data/Tiputini/trail_polylines.shp')
> plot(tbs)

So once you have loaded in a Shapefile, it is important to first find out what coordinate system is used, if any. This is done easily with the proj4string commands and functions in the package sp.

> tbs@proj4string  # This command shows us that there is no Coordinate System assigned to the TBS shapefile. 
## CRS arguments: NA

To be sure there are no coordinate data related to the file in any way, run the following code also from the sp package. Even though it seems there is no chosen projection, in this case our new projection will not be used correctly without using this command first.

> proj4string(tbs) <- NA_character_  # remove Coordinate Reference System information from TBS shapefile.

Now we can specify our new coordinate system, which will be the Projected Coordinate System UTM (Universal Transverse Mercator). This projects the earth on a flat surface and then uses x and y coordinates to pinpoint locations on a 2D surface. The command for this is CRS() from the sp package.

> proj4string(tbs) <- CRS("+proj=utm +zone=18 +south +datum=WGS84")  #Setting UTM WGS84 as our Coordinate Reference System (CRS).
> plot(tbs)

Importing and projecting .csv data

Now that we have defined the coordinate system of the shapefile, we can pull in the ranging data collected from the titi monkeys between July 2014 to June 2015. These data have been cleaned in advance and stored in a .csv file on GitHub.

> f <- curl("https://raw.githubusercontent.com/Callicebus/vignette/master/GPScoordinates.csv")
> points <- read.csv(f)
> head(points)
##   Obs.Sample.ID     Date Month        Group     UTMX    UTMY
## 1       OS92719 4-Jul-14     7 Callicebus L 371667.0 9929338
## 2       OS92719 4-Jul-14     7 Callicebus L 371647.7 9929378
## 3       OS92719 4-Jul-14     7 Callicebus B 371664.1 9929523
## 4       OS92719 4-Jul-14     7 Callicebus B 371747.5 9929529
## 5       OS92719 4-Jul-14     7 Callicebus B 371711.8 9929542
## 6       OS92719 4-Jul-14     7 Callicebus B 371700.1 9929542

As you can see, the coordinates in this file are already in the UTM projection, but we will have to tell R that. This is also done through commands in the sp package. First, you’ll have to convert the .csv file to a dataframe for spatial points, and then set the right coordinate system.

> spdf <- SpatialPointsDataFrame(coords = c(points[5], points[6]), data = points[4])  #Column 5 has the X-coordinates, column 6 has the Y-coordinates, and column 4 has important attribute data.
> str(spdf)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  1719 obs. of  1 variable:
##   .. ..$ Group: Factor w/ 3 levels "Callicebus B",..: 3 3 1 1 1 1 1 1 2 2 ...
##   ..@ coords.nrs : num(0) 
##   ..@ coords     : num [1:1719, 1:2] 371667 371648 371664 371747 371712 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "UTMX" "UTMY"
##   ..@ bbox       : num [1:2, 1:2] 371376 9929241 372130 9929694
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "UTMX" "UTMY"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr NA
> plot(spdf, pch = 1, col = spdf$Group)

> utm <- SpatialPoints(spdf, proj4string = CRS("+proj=utm +zone=18 +south +datum=WGS84"))
> 
> # If you would like to transform your data to Lat/Long rather then UTM, you
> # can use this code: latlong <- spTransform(utm, CRS('+proj=longlat
> # +datum=WGS84'))
> 
> utm <- SpatialPointsDataFrame(utm, data = points[4])
> plot(tbs)
> plot(utm, pch = 20, col = utm$Group, add = T)  #Simply adding add = TRUE to your code will overlay the second plot on the first.

If you want to zoom in or out on your plots, you can use a simple command from the zoom package. Here we have not put it in a code block, because you can only use this commend directly in your R console. Type the following:

inout.zoom()

Now you can zoom in by left-clicking on your plot. By left-clicking near the edge of your plot you can zoom out. Last, you can leave this function by hitting Escape.

NOTE: Always make sure that every file you work with has the same coordinate system. If you have different projections or datums, your layers will not overlay and you’ll likely feel like the baboon in this GIF. These problems are created by the fact that geographical and projected coordinate systems project data in a different way (in 3D versus 2D representations, respectively) and the fact that the Earth isn’t actually perfectly round!

Callicebus discolor

Part 2: Analyzing home range data

Creating home range polygons and estimating home range size

Now that we have all points plotted, we can draw polygons around the points of each group, representing the home ranges of the monkeys during the one-year period. We’ll do this with the help of a package build especially for home range measurements: adehabitatHR

NOTE: the functions in this package work best if you use the UTM projections!

First we calculate the utilization distribution with the help of the kernelUD() command. The ulitization distribution is a two dimensional probability density function that represents the probability of finding an animal in a defined area within its home range. Next we will use the command getverticeshr() to draw polygons around the kernel UD’s that have just been calculates.

> kernel <- kernelUD(utm[, 1], h = "href")
> TitiRange <- getverticeshr(kernel, 95)
> plot(TitiRange, border = TitiRange@data$id)
> plot(tbs, add = T)

There are many ways to calculate home range sizes, but by far the quickest and easiest way to calculate home range size is by using the as.data.frame() command. If you are working with Lat/Long coordinates, you will want to reproject them into UTM coordinates, because the function below only gives you a home range size in hectares if you use UTM coordinates (which are usually in meters).

> as.data.frame(TitiRange)
##                        id     area
## Callicebus B Callicebus B 6.195812
## Callicebus K Callicebus K 4.276863
## Callicebus L Callicebus L 6.643743

How are home ranges used?

Let’s look only at the home range of Callicebus group L. Using our points, we can also find out which parts of their home range are used more often, and which parts are used less frequently. First, let’s subset the data so we only have the ranging points of Callicebus L.

> library(RColorBrewer)
> Lpoints <- points[points$Group == "Callicebus L", ]

There are many ways to create a heat map, but the easiest way seems to be through ggplot. You use the ranging points from group Callicebus L and then create the heat map by adding stat_density2d(). You can add fancy colors with the RColorBrewer package.

> hothothot <- ggplot(Lpoints, aes(x = UTMX, y = UTMY)) + stat_density2d(aes(fill = ..density..), 
+     geom = "tile", contour = FALSE) + geom_point(size = 0.25) + scale_fill_gradientn(colours = rev(brewer.pal(7, 
+     "Spectral"))) + coord_fixed() + ggtitle("Titis Like it HOTHOTHOT")
> hothothot